%% FBMC/OQAM
%% FBMC-AWGN: BER verus SNR for FBMC/OQAM using polyphase structure

%%------------------------------------------------------------------------------------------------------------%%
%
%             "Prototype filter and filter bank structure"
% FBMC/OQAM with noise
%%------------------------------------------------------------------------------------------------------------%%
clc
clear all;
%%------------------------------------------------------------------------------------------------------------%%
%                                               Intialization
%%------------------------------------------------------------------------------------------------------------%%
nd=2; %number of OQAM sample per symbols (not change)
M=512;% number of sub-channels
K=4; % overlapping factor(good K=4)
lp=K*M-1;% prototype filter length
SNR= 0:5; % signal to noise ratio
nloop=10;  % Number of simulation loops
Rate0= zeros(1, length(SNR));
%%------------------------------------------------------------------------------------------------------------%%
%                                           Design of Prototype Filter
%%------------------------------------------------------------------------------------------------------------%%
%  low pass filter FIR
y=[1 0.97195983 sqrt(2)/2 0.23514695 ]; % coefficient of filter when K=4
yy=[1 0.911438 0.411438];% coefficient of filter when K=3
yyy=[1 sqrt(2)/2];% coefficient of filter when K=2
u=lp-1;
s=2*pi/(K*M);
for m=0:u-1
    r=m+1;
    p(m+1)=y(1,1)-2*y(1,2)*cos(r*s)+2*y(1,3)*cos(2*r*s)-2*y(1,4)*cos(3*r*s);% prototype filter equation K=4
    %p(m+1)=yy(1,1)-2*yy(1,2)*cos(r*s)+2*yy(1,3)*cos(2*r*s); % prototype filter equation K=3
    %p(m+1)=yyy(1,1)-2*yyy(1,2)*cos(r*s); % prototype filter equation K=2
    %p(m+1)=1; % prototype filter equation for OFDM
end
p=[0 p 0]; % extra delay sample (z^-2)
ty=p; %
%%------------------------------------------------------------------------------------------------------------%
%%------------------------------------------------------------------------------------------------------------%%
%%------------------------------------------------------------------------------------------------------------%%
%                                               Transmitter
%%------------------------------------------------------------------------------------------------------------%%
for ebn0=1:length(SNR)
    ebn0
    for iii=1:nloop
        % Data generation
        q=2;% binary data
        m=4;% modulation level
        In_Data=randi([0 q-1],M,1);% random input data M x 1
        % OQAM Modulation
        [OQAM_In_Data]=OQAM_preprocessing(In_Data,m,M); % OQAM modulation M x nd
        %%------------------------------------------------------------------------------------------------------------%%
        %                                               Transform Block
        %%------------------------------------------------------------------------------------------------------------%%
        %% beta carrier
        for k=0:M-1
            B(k+1)=(-1)^k*exp(-1*i*2*pi*k*(lp-1)/(2*M));
        end
        B=reshape(B,M,1);
        % ch1=Bi*Di   beta*data
        for k=0:M-1
            %   CHT1(k+1,:)=B(k+1)*Train_Data(k+1,:);
            CHD1(k+1,:)=B(k+1)*OQAM_In_Data(k+1,:);
        end
        %% 'ifft'
        CHD2=ifft(CHD1);
        %%------------------------------------------------------------------------------------------------------------%%
        %                                              Polyphase Filters SFB
        %%------------------------------------------------------------------------------------------------------------%%
        % polyphase filters SFB
        for k=0:M-1
            j(k+1,:)=p(k+1+M*(0:K-1));%% Aq coefficients
            ss=upsample(j(k+1,:),2);
            CHD3(k+1,:)=cconv(CHD2(k+1,:),ss,nd);%% ch2*Aq
        end
        %%------------------------------------------------------------------------------------------------------------%%
        %                                              P/S Conversion
        %%------------------------------------------------------------------------------------------------------------%%
        % Up sample and delay chain
        % instead of up sample and delay chain
        [a b]=size(CHD3);
        CHD4=reshape(CHD3,a*b,1);
        
        %%------------------------------------------------------------------------------------------------------------%%
        %                                              Channel Effect
        %%------------------------------------------------------------------------------------------------------------%%
        %% channel effects
        CHD5=reshape(CHD4,a*b,1);
        %%
        CHD6=CHD5;
        % CHD7 = awgn(CHD6,ebn0,'measured');
        signal_in_dB=10*log10(std(CHD6)^2);
        noise_in_dB=signal_in_dB-SNR(ebn0);
        noise=(10^(noise_in_dB/10)^(1/2))*randn(size( CHD6, 1),size( CHD6, 2));
        CHD7=CHD6+noise;
        
        %%------------------------------------------------------------------------------------------------------------%%
        %                                              Receiver
        %%------------------------------------------------------------------------------------------------------------%%
        %%------------------------------------------------------------------------------------------------------------%%
        %                                             S/P Conversion
        %%------------------------------------------------------------------------------------------------------------%%
        CHD7=reshape(CHD7,a,b);
        %%------------------------------------------------------------------------------------------------------------%%
        %                                              Polyphase Filters AFB
        %%------------------------------------------------------------------------------------------------------------%%
        for k=0:M-1
            ee=upsample(j(M-k,:),2);
            CHD8(k+1,:)=cconv(CHD7(k+1,:),ee,nd);
        end
        %%------------------------------------------------------------------------------------------------------------%%
        %                                             Transform Block
        %%------------------------------------------------------------------------------------------------------------%%
        %  'fft'
        CHD9=fft(CHD8);
        % Beta demodulator
        for k=0:M-1
            BB(k+1)=(-1)^k*exp(1*i*2*pi*k*(lp-1)/(2*M));
        end
        BB=reshape(BB,M,1);
        %
        for k=0:M-1
            CHD10(k+1,:)=BB(k+1)*CHD9(k+1,:);
        end
        %%------------------------------------------------------------------------------------------------------------%%
        %
        %%------------------------------------------------------------------------------------------------------------%%
        
        %%------------------------------------------------------------------------------------------------------------%%
        %OQAM Demodulation
        [Out_Data]=OQAM_postprocessing(reshape(CHD10,M,nd),m,M);
        
        %%------------------------------------------------------------------------------------------------------------%%
        %                                              BER
        %%------------------------------------------------------------------------------------------------------------%%
        
        [err0, rate0]= symerr( Out_Data,In_Data);
        
        Rate0(ebn0)= Rate0(ebn0) + rate0;
        
    end
    % Average value
    
    Rate0(ebn0)= Rate0(ebn0)/nloop; %
    
end
%%------------------------------------------------------------------------------------------------------------%%
%                                              Figure
%%------------------------------------------------------------------------------------------------------------%
f1 = figure(1);
set(f1,'color',[1 1 1]);
hold on
semilogy(SNR,Rate0,'r-^')
xlabel( 'Signal-to-Noise Ratio(SNR) in [dB]')
ylabel( 'Bit Error Rate (BER)')
